Dynamic interspecies interactions and robustness in a four‐species model biofilm

Abstract Interspecific interactions within biofilms determine relative species abundance, growth dynamics, community resilience, and success or failure of invasion by an extraneous organism. However, deciphering interspecific interactions and assessing their contribution to biofilm properties and function remain a challenge. Here, we describe the constitution of a model biofilm composed of four bacterial species belonging to four different genera (Rhodocyclus sp., Pseudomonas fluorescens, Kocuria varians, and Bacillus cereus), derived from a biofilm isolated from an industrial milk pasteurization unit. We demonstrate that the growth dynamics and equilibrium composition of this biofilm are highly reproducible. Based on its equilibrium composition, we show that the establishment of this four‐species biofilm is highly robust against initial, transient perturbations but less so towards continuous perturbations. By comparing biofilms formed from different numbers and combinations of the constituent species and by fitting a growth model to the experimental data, we reveal a network of dynamic, positive, and negative interactions that determine the final composition of the biofilm. Furthermore, we reveal that the molecular determinant of one negative interaction is the thiocillin I synthesized by the B. cereus strain, and demonstrate its importance for species distribution and its impact on robustness by mutational analysis of the biofilm ecosystem.

In contrast to typical laboratory conditions of growth in liquid culture, bacteria in natural environments and those contaminating hospitals, or industrial and food-processing procedures are more often found in multicellular surface-associated communities known as biofilms (Costerton et al., 1987;Flemming et al., 2016;Hall-Stoodley et al., 2004). Such biofilms are generally complex communities harboring numerous bacterial species in close spatial proximity (Elias & Banin, 2012). Diverse physical and social interactions between species take place in these communities. They are considered to determine not only the structure and spatial organization of the biofilm but also its global functions by modulating gene expression in the different species (Bridier et al., 2017;Burmølle et al., 2014;Liu et al., 2016;Rendueles & Ghigo, 2012). The physiology of each microbial species in complex, multispecies biofilms might be distinct from that in monospecific biofilms (L. B. S. Hansen et al., 2017;Liu et al., 2019). Moreover, multispecies biofilms exhibit emergent properties such as increased tolerance against antimicrobial agents (Bridier et al., 2012;Burmølle et al., 2006;Schwering et al., 2013;Yan & Bassler, 2019), synergistic degradation of toxic compounds (Breugelmans et al., 2008;Perera et al., 2019;Yoshida et al., 2009), stronger defense against protozoan grazing (Koh et al., 2012;Raghupathi et al., 2018), increased virulence in infection (Pastar et al., 2013;Wang et al., 2020) and protection against the action of biocides Sanchez-Vizuete, Orgaz, et al., 2015;Yan & Bassler, 2019). Studies on multispecies biofilms have also reported enhanced stress resistance, productivity, or biomass production (Burmølle et al., 2006;Lee et al., 2014;Liu et al., 2019;Ren et al., 2015), and, importantly, "community-intrinsic properties" (Madsen et al., 2018) emerging from the social interactions between members of the biofilm and which may be important for its interaction with its environment.
It is thus important to decipher these interactions, positive or negative, at the molecular and biochemical levels to better understand the ecological and evolutionary factors that drive community function in natural or engineered systems Ziesack et al., 2019). However, the number and types of interactions within multispecies biofilms are expected to grow very rapidly with the number of species present in the biofilm . Characterization of the interactions in complex biofilms and their underlying molecular mechanisms remains a challenge, as well as the evaluation of the importance of these interactions for the overall robustness of the structure and functionalities of these biofilms (Røder et al., 2020).
Here, we constructed a biofilm community from four species isolated from a biofilm consortium contaminating a milk processing plant, to analyze interspecies interactions and robustness to environmental stresses. This multispecies biofilm was highly reproducible, allowing us to use it as a model to study the dynamic interactions that take place between the species during its development. We were thus able to test the resistance of this complex biofilm and its formation process towards different continuous or transient perturbations. Then we identified the molecular determinant and assessed the contribution of one major interspecies interaction to the overall robustness of the biofilm.

| Bacterial strains and culture conditions
Bacillus cereus ATCC14579 T and B. cereus ATCC10987 were obtained from American Type Culture Collection (ATCC; https://www. lgcstandards-atcc.org). The 13 other bacterial strains used in this study (following paragraph and Table 1) were isolated from a biofilm formed on a gasket in a milk pasteurization line (Mettler & Carpentier, 1997

| Experimental model biofilm formation
The biofilm development model of Maris (Maris, 1992)  Goodfellow). Before use, they were prepared as previously described (Leriche & Carpentier, 1995). Cultures in the mid-exponential growth phase were diluted into a fresh TSB medium to give a mixture con-

| Planktonic coculture
Strains were individually precultivated in TSB at 30°C with agitation, harvested by centrifugation, and resuspended in 1/20 TSB.
They were diluted and mixed as above to give 5 × 10 6 cfu ml −1 of each strain in 20 ml 1/20 TSB. This planktonic coculture mixture was placed in a 9-cm diameter Petri dish and incubated for 24, 48, or 72 h at 30°C with agitation at 60 rpm to prevent biofilm formation.

| Electron microscopy
Scanning electron microscopy (SEM) experiments were conducted as described (Couvigny et al., 2017). Images were acquired and analyzed at the Microscopie et Imagerie des Micro-organismes, Animaux et Aliments (MIMA2) microscopy and imaging platform.

| Estimation of growth parameters with a logistic model
To represent the growth of a single-species population, we used a simple logistic equation. Adding a latency phase with a constant lag time = τ, the differential equation reads: dynamics, with population x t ( ) converging towards K (the carrying capacity). Integration of this differential equation leads to the wellknown logistic function: Given an experimental growth curve T X ( , ) i i , we set the initial condition x 0 to be the bacterial count immediately after the 90 min adhesion step and we estimated parameters K, μ, and τ to minimize the following (constrained) least-square error:

| Statistical analysis
Values for end-point biofilm growth experiments are expressed as the mean ± standard deviation calculated from at least five independent experiments. For statistical analysis, cfu for each strain was expressed as logarithms, except where noted. Statistical significance of measured differences was determined using a two-way analysis of variance with Tukey's posttest as implemented in "R". For Figure 6, in which members were omitted from the biofilm community, the significance of changes in the growth of the species was estimated by the Mann-Whitney Wilcoxon test, compared to growth in a four-species biofilm. In the particular case of biofilms growing in non-diluted medium (Figures 4 and 8), cfu for Rhodocyclus were in all cases below the limit of detection, comparison with biofilms grown in other conditions was made by assigning random values between 10 2 and 10 4 cfu cm −2 (i.e., the detection limit plus or minus 1 on a log scale) to the Rhodocyclus measurements for each of the seven data sets obtained from biofilms grown in a non-diluted medium. The procedure was repeated twenty times and the p values retained are the upper 95% confidence levels. No significant difference resulted in the estimated significance values when the range of random values assigned to the Rhodocyclus cfu was increased tenfold, nor when the values were reduced further.

| Construction of a multispecies, model biofilm based on a natural ecosystem
To construct a model, multispecies biofilm we exploited an industrial biofilm consortium, isolated from an industrial food preparation device (Mettler & Carpentier, 1997), comprising 13 strains corresponding to 7 species belonging to 5 genera (Table 1). Based on their capacity to form monospecies biofilms under laboratory conditions, we chose the best biofilm formers from each genus to simplify the consortia and facilitate subsequent analysis (see Section 2 and Table 1). Among these selected strains, we substituted the Bacillus strain with a wellstudied and genetically tractable laboratory strain (B. cereus ATCC14579 T ). We then tested combinations of the strains, including a representative of each genus, for their ability to form a multispecies biofilm. None of the Staphylococcus strains persisted at measurable levels in the laboratory multispecies biofilms.
Hence, the consortium selected for further analyses contained the following strains: Rhodocyclus sp. CCL5, P. fluorescens CCL49, K. varians CCL56, and B. cereus ATCC14579 T (henceforth referred to, for simplicity, as Rhodocyclus, P. fluorescens, K. varians, and B. cereus, respectively) thus representing four of the five genera initially present in the industrial biofilm.

| Community structure in biofilm versus planktonic lifestyles
To determine whether the steady state reached by the four-species community is biofilm-specific, we compared the population distributions in biofilm and liquid culture ( Figure 3). The biofilm population structure (Figure 3a), with a majority of Rhodocyclus (73.6%) and P. fluorescens (25.2%) and considerably less of K. varians (0.2%) and B. cereus (1.0%) was markedly different from that of the planktonic grown community ( Figure 3b). Proportions in liquid medium were more even, with B. cereus, Rhodocyclus, P. fluorescens, and K.
varians comprising respectively 38.2%, 29.8%, 18.0%, and 14.0% of the total cell counts. Since surface-associated biofilm populations release cells to the bulk liquid, and also recruit cells from the planktonic phase (Houry et al., 2012) it was interesting to determine whether the bacterial community composition present in the liquid phase above the biofilm was characterized by a biofilm-like or by a planktonic-like species distribution. As seen in Figure 3c, the composition of the biofilm supernatant community (Rhodocyclus, 62.6%; P. fluorescens, 34.6%; B. cereus, 2.8%; and K. varians, 0.01%) was more similar to that of the biofilm than to that observed in the pure planktonic culture. This suggests that, in the experimental setup, exchanges between biofilm and supernatant are more influenced by the seeding of the planktonic phase from the biofilm during the 24 h between washes than by the migration of cells from the planktonic phase to the biofilm. That is, the biofilm species composition is a F I G U R E 1 Evolution of species abundances during biofilm development. Times are in hours, after the 90-min adhesion phase. Results are averages of at least five independent experiments; error bars correspond to one standard deviation. The evolution of the biofilm composition is analyzed in terms of changes of each species' abundance compared to that of the previous reading: *p < .05, **p < .01, ***p < .01 (Tukey's HSD). Bc, B. cereus; HSD, honestly significant difference; Kv, K. varians; Pf, P. fluorescens; Rh, Rhodocyclus F I G U R E 2 Physical Interactions between species in the biofilm. Scanning electron micrograph of the four-species biofilm at T = 72 h (a) and of a biofilm grown under identical conditions but containing a Bacillus cereus ΔtclE-H mutant lacking thiocillin production in place of the wild-type strain (b). The four species are indicated by arrows in panels a and b: (1) Rhodocyclus, (2) Pseudomonas fluorescens, (3) Kocuria varians, (4) B. cereus. Note that in biofilms containing the wild-type B. cereus, K. varians numbers are much reduced, such that none are observed in the image (a) which is representative of the majority of the observed microscopic fields. Scales are shown at the bottom right property of the biofilm and not a reflection of the equilibrium attained in the planktonic phase. Thus, the equilibrium reached in biofilm appears specific to this mode of growth, and these results support the hypothesis that interspecies interactions are different in the biofilm compared to the planktonic culture, driving the system to a different equilibrium population distribution.

| Robustness of the multispecies biofilm to perturbations
To evaluate the robustness of the biofilm, we applied transient or continuous changes of experimental conditions, and then compared the relative species abundances and the total cell numbers of the biofilm with those obtained under control standard conditions (see Section 2).
Neither total population counts nor species abundances at steady state were significantly altered by initial transient perturbations modifying either the global physiological state or the composition of the starting inoculum ( Figure 4). In particular, drastic reduction in the inoculum (from 10 5 to 10 2 cfu ml −1 ) of B. cereus, the laboratory strain that has been substituted to the original Bacillus strain, did not change the relative species proportions and B. cereus was in all cases able to attain equilibrium cell densities (∼6.10 5 cfu cm −2 ) demonstrating its capacity to grow in the biofilm. In contrast, two continuous perturbations did change the final structure of the community: Changes in the substrate surface properties (glass in place of steel) affected the equilibrium proportion of B. cereus, whose proportion increased from 1.0% to 16.3%, while the 20-fold concentration of the culture medium affected both the biomass and species distribution of the biofilm (Figure 4). With this latter perturbation, total cell numbers increased threefold and relative species abundances were significantly modified. While P. fluorescens reached 93.3% of the total bacterial counts at 6.2 × 10 7 cfu cm −2 , B. cereus cell counts were increased 18-fold to 4.3 × 10 6 cfu cm −2 (6.5% of the total cell counts), K. varians remained a minority at 1.8 × 10 5 cfu cm −2 (0.3%) and Rhodocyclus, the major component of the community under standard conditions, was no longer detectable (less than 1000 cfu cm −2 ). Less drastic modification of the culture medium (twofold concentration) had no measurable effect.
These results show that the establishment of the multispecies biofilm is resistant to initial, even quite considerable, transient To gain insights into the mechanisms underlying these effects, we fitted a simple growth model to the experimental data and estimated growth parameters for the different species in single-and four-species biofilms (Table 2). We used least-square minimization to estimate three parameters: the maximal growth rate μ, the carrying capacity K and the lag time τ (see Section 2). For Rhodocyclus in a single-species biofilm, the fit was not possible due to initial counts below the limits of detection ( Figure 5a); while for B. cereus it was not relevant as the net growth rate within the biofilm is zero, both alone and in the biofilm community (Figure 5d). In the case of K. varians, the considerably reduced carrying capacity when present as part of a four-species biofilm, as compared to its growth a single-species biofilm, and its apparent lack of growth in the mixed community suggests that one or more of the other species have a negative influence on its growth. In the case of P. fluorescens, a negative effect is also observed, although of a different nature: parameter estimates indicated that growth rate and carrying capacity for P. fluorescens were comparable in the two cultivation conditions, whereas its lag time increased greatly when it was cultivated in the four-species biofilm ( Table 2). Since the adhesion step resulted in similar P. fluorescens populations at time zero (Figure 5b), it is apparent that P. fluorescens development was specifically delayed in the fourspecies biofilm.
In the case of Rhodocyclus, differences were seen in the initial development stages (Figure 5a). Adherent bacteria are undetectable in single-species conditions (less than 1000 per square centimeter, compared with 3000 in the four-species mixture), and remain undetectable after 24 h of growth. However, the population densities after 48 or 72 h are indistinguishable in the single-species and fourspecies mature biofilms. Thus, the presence of the three other species facilitates adhesion and initial development but does not increase the carrying capacity for Rhodocyclus.
Population densities of B. cereus at each biofilm development stage were similar in monospecies and four-species biofilms (Figure 5d).

F I G U R E 4
Effects of experimental perturbations on species distribution in the biofilm. Stationary phase: Inoculation with stationary-phase cultures instead of exponential-phase cultures; Bc 1/10, Bc 1/100, or Bc 1/1000: reduction of the amount of B. cereus in the initial inoculum from 5 × 10 6 to 5 × 10 5 , to 5 × 10 4 , or 5 × 10 3 cfu ml −1 , respectively; glass coupon: replacement of stainless steel chip with glass; non-diluted growth medium: growth in undiluted TSB instead of 1/20 TSB. Results are averages of five or more independent experiments; error bars correspond to one standard deviation. *p < .05, ***p < .001, significance (Tukey's HSD) of differences in species proportions relative to those in biofilms grown under standard conditions. Bc, B. cereus; HSD, honestly significant difference; Kv, K. varians; Pf, P. fluorescens; Rh, Rhodocyclus; TSB, tryptic soy broth Though the net growth rate of B. cereus in the biofilm under standard conditions is zero, this species is capable of multiplication within the biofilm since, when it is inoculated at a level of 2.5 × 10 3 cfu cm −2 , the final density nevertheless reaches 10 5 (data for "Bc 1/1000" in Figure 4).
Hence, B. cereus is neither favored nor disadvantaged in the multispecies biofilm.
These results reveal the existence of distinct types of interaction, both negative and positive between members of the multispecies F I G U R E 5 Bacterial growth in single-species and four-species biofilms. Cfu measured at 24-h intervals in one-species (dotted lines) and four-species (solid lines) biofilms for (a) Rhodocyclus, (b) Pseudomonas fluorescens, (c) Kocuria varians, and (d) Bacillus cereus. Error bars correspond to one standard deviation; results are the averages of at least five independent experiments. cfu for Rhodocyclus growing as a single-species biofilm were below the limits of detection at 0 and 24 h T A B L E 2 Estimation of K. varians and P. fluorescens growth parameters in single-and four-species biofilm using a logistic growth model

Single-species biofilm
Four-species biofilm Rhodocyclus sp. Kocuria varians 1.93 2.1 × 10 7 3 -1.0 × 10 5 -Bacillus cereus -7.6 × 10 4 --2.3 × 10 5 -Note: The symbol "-" indicates that the estimation of µ and was inconclusive on the corresponding data set, either because of insufficient information at early time points (Rhodocyclus in single-species biofilm) or because of an absence of apparent growth between t = 0 and 72 h (K. varians in four-species biofilm and B. cereus in single-and four-species biofilms), see Figure 4d. In those cases, the value indicated for the carrying capacity K is the mean value of the population in the corresponding data set.
biofilm. To further investigate these interactions, we assessed the effects of species omissions and substitutions on the final composition of the multispecies biofilm.

| Effects of species omissions on biofilm composition
We 3.7 | Identification of the B. cereus product responsible for its negative interaction with K. varians

Colonies of B. cereus induced clear zones of inhibition of K. varians
when these two strains were grown together on nutrient agar plates (Figure 7a), suggesting that this strain produces one or more inhibitory substances. Comparing the activities of different B. cereus strains, we noted that the strain ATCC10987 does not inhibit K.
varians growth. Comparison of its published genome with that of the laboratory strain used in the present study, B. cereus ATCC14579, indicated a group of genes BC5071 to BC5102 which, except for one, were absent from ATCC10987. Further sequence comparisons confirmed that this region is highly variable from strain to strain in the B. cereus-B. thuringiensis-B. anthracis group ( Figure A1). This region is part of the defined biosynthetic gene cluster tclA-X (BC5094 to BC5071) which is involved in the production of thiocillins, modified peptide antibiotics Wieland Brown et al., 2009).
Deletion of the four structural genes coding the thiocillin precursor,

| Molecular mechanism of interference competition by growth inhibition
We found that the negative effect on K. varians in the four-species biofilm was due to the production of thiocillin by B. cereus and that this negative interaction was partly mitigated by the presence of P. fluorescens and Rhodocyclus. This bacteriocin is active against Gram-positive bacteria, and B. cereus was found to inhibit the growth of all Staphylococcus tested, explaining the failure of these strains to form a mixed-species biofilm ( Table 1) The coexistence of antagonistic strains in stable, or in cyclically evolving, proportions is predicted by theory for a range of interaction parameters (Chesson, 2000;Czárán et al., 2002;Hassell et al., 1994) and has been demonstrated in laboratory experiments (Kerr et al., 2002). Notably, these latter authors showed that antagonistic strains can coexist if interactions and dispersal occur at a local scale, coexistence under these conditions being associated with structured spatial organization of the strains. This biofilm structuration, which can protect sensitive species from the antagonistic effects of others (Kim et al., 2011), is also favored by positive interactions (Breugelmans et al., 2008;S. K. Hansen et al., 2007), such as the observed stimulation of Rhodocyclus adhesion by the biofilm community ( Figure 5d). Nevertheless, the mechanism of mitigation by P. fluorescens and Rhodocyclus of the inhibitory action of B. cereus on K. varians may involve nonspecific protection by the physical presence of these two species-for example, physical separation and/or adsorption of the thiocillin molecule at the bacterial surface.

| Interspecies interactions and robustness of the biofilm community
The establishment of the four-species biofilm was resistant in the face of the transient perturbations that we tested, but the species composition changed in response to continuous perturbations, where the environment was permanently modified. We then looked more closely at the effect on biofilm community stability of one of the interspecies interactions, in particular, the negative effect of B. cereus on the growth of K. varians. As described above, the substitution of the B. cereus wild-type by the thiocillin mutant strain altered the final composition of the biofilm, permitting increased growth of K. varians which became codominant together with Rhodocyclus and P. fluorescens. Neither the wild-type nor the mutant biofilm was robust to strong continuous perturbations (compare Figures 4 with 8). Moreover, the mutant biofilm appeared to be less resistant than the wildtype to transient perturbations, suggesting that in this model the major negative interaction may play a role in robustness that is only partially compensated by other stabilizing interactions (Burmølle et al., 2006;Lee et al., 2014). This is in agreement with the results of a study by Thompson et al. (2020) concerning a community composed of bacteria isolated from a potable water distribution system, where the authors detected redundant interspecies interaction effects. The interaction between one species and either one of two others had a positive effect on the biomass of their model biofilm, and this effect was more marked in the absence of a fourth, whose presence independently compensated for the loss of the interactions. Redundant interactions were also brought to light in our model biofilm but it must be remembered that the bacteria in the biofilm have had little time to coevolve together-probably 15-20 generations-and that the major negative interaction, though it may be important in this model, is not part of an evolved ecosystem.
Growth of the biofilm in the undiluted medium resulted in considerably higher numbers of B. cereus (Figure 4), possibly due to a capacity for exploiting the increased nutrient availability, to changed interactions in the biofilm community, or a combination of the two effects. In contrast, cell density of the B. cereus ΔtclE-H was decreased in rich medium (Figure 8), both in relation to growth under standard conditions and to the wild-type in the undiluted medium.
These observations suggest that the differential reaction of the wildtype and mutant biofilms to growth in the undiluted medium is at least partly explained by a modification of interspecies interactions other than competition.

| CONCLUSIONS
Emergent properties of bacterial communities grown as biofilms, driven by social interactions, have huge implications for research and practical knowledge in such contexts as human health, food safety, rhizosphere role in plant growth, or even bioremediation.
One approach to understanding these social interactions is to create and study artificial biofilm consortia in the laboratory.
However, very few studies report such reconstructions of multispecies biofilm and elucidate the interspecies interaction networks that take place within. Moreover, the molecular determinant of these interactions and the analysis of their impact on the biofilm ecosystem properties have been reported in only a few studies.
Here, we not only deciphered the active network of interactions that shapes a four-species biofilm community and determine its robustness but also identified the molecular determinant of one of these interactions and revealed how it impacts the structure and properties of this community.

APPENDIX
F I G U R E A1 Alignments of genomic regions around the locus containing the thiocillin biosynthesis genes in Bacillus cereus ATCC14579. Genomes aligned are the two B. cereus strains used in this study ATCC14579 T and ATCC10987, and the representative genomes as defined at NCBI of other B. cereus group clades. Darker gray trapezoids join sequence blocks with over 90% base identity between strains; the lighter gray regions are at the position of a large and variable CDS (BC5055 in strain ATCC14579) containing peptidase and multiple mucin-binding domains. CDS are colored light or dark yellow to distinguish transcriptional units as described by (Kristoffersen et al., 2012), Genome Biology, 13:R30); outside the annotated thiocillin production genes (red and pink), white indicates that no transcripts were detected in the latter study. Notable groups of genes in B. cereus ATCC14579 are marked: (a) thiocillin structural genes, four identical short ORFs, BC5090 to BC5087, (b) lantibiotic synthesis enzyme genes, BC5086 to BC5083, (c) polyketide resistance, BC5091 and BC5092, (d) macrolide export, BC5072 and BC5071, (e) cytolysin precursor BC5101 (present in a number of B. cereus group strains, elsewhere on the chromosome of B. anthracis Ames, 70% predicted amino acid identity to perfringolysin O) and (f), outside the island, siderophore binding and import, BC5106 to BC5103. The island in B. anthracis Ames is an integrated prophage